13_Pooled_Results

## Load functions to create posterior distributions 

# to generate posteriors distributions from beta distributions using
# noninformative priors within tibbles

source(here("final_analyses", "script", "functions", # file path
            "noninformative_posterior.R"))           # file name

# To generate normal conjugate posterior probabilities

source(here("final_analyses", "script", "functions",           # file path
            "conjugate_normal_posterior.R"))                   # file name

# To generate posterior probabilities with objects from rma()
# (it requires "conjugate_normal_posterior.R" )

source(here("final_analyses", "script", "functions", "pooled", # file path
            "pooled_normal_approximation.R"))                  # file name

# To generate posterior probabilities with multiple priors
# (it requires "conjugate_normal_posterior.R" )

source(here("final_analyses", "script", "functions", "pooled", # file path
            "pooled_normal_approximation_multiple_priors.R"))  # file name

# To compile all posteriors into a tibble

source(here("final_analyses", "script", "functions", "pooled", # file path
            "pooled_tibble_all_posteriors.R"))                 # file name

## Load functions to plot

# To create summary table

source(here("final_analyses", "script", "functions", "pooled", # file path
            "pooled_summary_table_death.R"))                   # file name

# To create summary table

source(here("final_analyses", "script", "functions", "pooled", # file path
            "pooled_summary_table_discharge.R"))                   # file name

# To plot risk difference distribution

source(here("final_analyses", "script", "functions", "pooled", # file path
            "pooled_risk_difference_plot.R"))                  # file name

# To plot cumulative probability of risk difference

source(here("final_analyses", "script", "functions", "pooled", # file path
            "pooled_cumulative_risk_difference_plot.R"))       # file name

# To plot prior, data and posterior distributions side-by-side

source(here("final_analyses", "script", "functions", "pooled", # file path
            "pooled_data_prior_posterior_plot.R"))             # file name

# To plot the posterior distribution

source(here("final_analyses", "script", "functions", "pooled", # file path
            "pooled_posterior_difference_plot.R"))             # file name

# To plot the cumulative probability of the posterior distribution

source(here("final_analyses", "script", "functions", "pooled",  # file path
            "pooled_cumulative_risk_difference_plot.R"))        # file name

# To plot the cumulative probability of the posterior distribution

source(here("final_analyses", "script", "functions", "pooled",  # file path
            "pooled_posterior_cumulative_plot.R"))              # file name

# Original data extraction table

d = readxl::read_excel(here("final_analyses", "data", # file path
                            "extracted-data.xlsx"),   # file name
                       na = "NA")   

# Change name of columns and remove columns related to symptoms
d = clean_names(d) %>% 
  select(-starts_with("symptom"))

Overview

During the data extraction process, we realized that some trials did not provide data on separate subgroups. Instead, they would only provide pooled data. Let’s check how much data we actually have.

How many distinct trials do we have regarding mortality?

d %>% 
  filter(outcome == "mortality") %>% 
  distinct(trial) %>% 
  flextable() %>% 
  autofit()

How many distinct trials do we have regarding hospital discharge?

d %>% 
  filter(outcome == "discharge") %>% 
  distinct(trial) %>% 
  flextable() %>% 
  autofit()

How much data do we have on “pooled subgroups”?

d %>%
  filter(str_detect(oxygen, "pooled"),
         outcome == "mortality") %>% 
  group_by(outcome, trial, oxygen) %>%
  arrange(oxygen) %>%
  count() %>%
  select(-n) %>% 
  rename("subgroup" = oxygen) %>% 
  flextable() %>%
  autofit()

d %>%
  filter(str_detect(oxygen, "pooled"),
         outcome == "discharge") %>% 
  group_by(outcome, trial, oxygen) %>%
  arrange(oxygen) %>%
  count() %>%
  select(-n)%>% 
  rename("subgroup" = oxygen) %>% 
  flextable() %>%
  autofit()

What distinct trials have data on pooled and separate subgroups regarding mortality?

d %>%
  mutate(pooled = oxygen) %>% 
  filter(str_detect(pooled, "pooled"),
         outcome == "mortality") %>% 
  select(trial, pooled)%>% 
  distinct(trial) %>% 
  rename("pooled" = trial)%>% 
  flextable() %>% 
  autofit()

d %>%
  filter(!str_detect(oxygen, "pooled"),
         outcome == "mortality") %>% 
  select(trial, oxygen) %>% 
  distinct(trial) %>% 
  rename("separate" = trial) %>% 
  flextable()%>% 
  autofit()

What distinct trials have data on pooled and separate subgroups regarding hospital discharge?

d %>%
  mutate(pooled = oxygen) %>% 
  filter(str_detect(pooled, "pooled"),
         outcome == "discharge") %>% 
  select(trial, pooled)%>% 
  distinct(trial) %>% 
  rename("pooled" = trial)%>% 
  flextable() %>% 
  autofit()

d %>%
  filter(!str_detect(oxygen, "pooled"),
         outcome == "discharge") %>% 
  select(trial, oxygen) %>% 
  distinct(trial) %>% 
  rename("separate" = trial) %>% 
  flextable()%>% 
  autofit()

Mortality outcome

Now, let’s pool all the data and analyze the mortality outcome. I believe odds ratio is the most appropriate effect size measure since there is great heterogeneity in baseline risk in each subgroup.

Prior

not_recovery = 
  d %>%
  filter(trial != "RECOVERY") %>% 
  mutate(oxygen = case_when(
    oxygen == "simple oxygen" ~ "simple_oxygen",
    oxygen == "non-invasive ventilation" ~ "non_invasive",
    oxygen == "invasive ventilation" ~ "invasive",
    oxygen == "pooled simple or non-invasive ventilation" ~ "POOLED_simple_non-invasive",
    oxygen == "pooled non-invasive or invasive ventilation" ~ "POOLED_non-invasive_invasive",
    oxygen == "pooled simple, non-invasive or invasive ventilation" ~ "POOLED_simple_non-invasive_invasive"
  )) %>% 
  unite(trial_subgroup, c(trial, oxygen))


not_recovery %>% 
  filter(outcome == "mortality") %>% 
  select(-outcome) %>% 
  flextable() %>% 
  autofit()

Version 1: Random-Effect Meta-Analysis

prior_death =
  rma(
    # Outcome
    measure = "OR", # Odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = not_recovery %>% filter(outcome == "mortality"),

    slab = paste(trial_subgroup),
    method = "REML"
  )

# Visualize it
forest(prior_death,
       showweights = T,
       atransf=exp)

### Using escalc()

# Same results as rma() with random-effects

# effect_sizes =
#   escalc(
#     # Outcome
#     measure = "OR", # risk difference
#     
#     # Tocilizumab group
#     ai = events_toci,
#     n1i = total_toci,
# 
#     # Control group
#     ci = events_control,
#     n2i = total_control,
#     
#     data = not_recovery %>% filter(outcome == "mortality"),
# 
#     slab = paste(trial_subgroup)
#   )
# 
# res <- rma(yi, vi, data=effect_sizes) 
# 
# forest(res,
#        atransf=exp)

Version 2: Summing up all events + Fixed-Effect MA

total_death = not_recovery %>% 
  filter(outcome == "mortality") %>% 
  # Sum all events
  mutate(across(events_toci:total_control, ~ sum(.x))) %>% 
  # As all rows are now the same, pick 1
  slice(1) %>% 
  # Rename trial as "Total"
  mutate("trial_subgroup" = "Total")

prior_death_v2 =
  rma(
    # Outcome
    measure = "OR", # Odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = total_death,

    slab = paste(trial_subgroup),
    method = "REML"
  )

# Visualize it
forest(prior_death_v2,
       atransf=exp)

RECOVERY

recovery = 
  d  %>% 
  filter(trial == "RECOVERY",
         oxygen != 'NA') %>% 
  mutate(oxygen = case_when(
    oxygen == "simple oxygen" ~ "simple_oxygen",
    oxygen == "non-invasive ventilation" ~ "non_invasive",
    oxygen == "invasive ventilation" ~ "invasive")) %>% 
  unite(trial_subgroup, c(trial, oxygen))

recovery %>% 
  filter(outcome == "mortality") %>% 
  select(-outcome) %>% 
  flextable() %>% 
  autofit()

Version 1: Random-Effect Meta-Analysis

recovery_death =
  rma(
    # Outcome
    measure = "OR", # Odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = recovery %>% filter(outcome == "mortality"),

    slab = paste(trial_subgroup),
    method = "REML"
  )

# Visualize it
forest(recovery_death,
       showweights = T, 
       atransf=exp)

Version 2: Summing up all events + Fixed-Effect MA

total_death = recovery %>% 
  filter(outcome == "mortality") %>% 
  # Sum all events
  mutate(across(events_toci:total_control, ~ sum(.x))) %>% 
  # As all rows are now the same, pick 1
  slice(1) %>% 
  # Rename trial as "Total"
  mutate("trial_subgroup" = "Total")

recovery_total_death =
  rma(
    # Outcome
    measure = "OR", # Odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = total_death,

    slab = paste(trial_subgroup),
    method = "REML"
  )

# Visualize it
forest(recovery_total_death,
       atransf=exp)

Posterior distribution

Data overview

Regarding the mortality outcome, odds ratio lower than 1 means benefit.

I am going to use data from Version 1 analyses. The effect size here will be log-odds ratio.

Regarding the skeptical prior, its mean is centered in 0, which is equal to log(1).

Regarding the optimistic prior, its mean is centered in -0.223, which is equal to log(0.8). I chose the value 0.8 because the RECOVERY trial was designed to detect a relative risk reduction of 20%.

Regarding the pessimistic prior, its mean is centered in 0.223, which is equal to log(1.25). I mirrored the optimistic prior’s mean in the log-odds ratio scale.

multiple_priors_death %>% flextable() %>% autofit()

all_posteriors_death = pooled_tibble_all_posteriors(
  
  # Output from rma() with only RECOVERY
  recovery_death, 
  
  # Output from normal_approximation_multiple_priors()
  multiple_priors_death 
)

Summary

# X axis
xlim_inf = 0.6
xlim_sup = 1.1
ticks_spacing = 0.1

# Assure subgroup order
# https://stackoverflow.com/questions/12774210/
# how-do-you-specifically-order-ggplot2-x-axis-instead-of-alphabetical-order

level_order = c("Pessimistic",
                "Optimistic",
                "Skeptical",
                "Evidence-based",
                "Non-informative")

p1 = all_posteriors_death %>%
  
  # make it tidy for ggplot
  pivot_longer(
    cols = c("Non-informative":"Pessimistic"),
    names_to = "dist", values_to = "samples"
  ) %>%
  
  # exp() to transform log-odds ratio into Odds Ratio
  
  mutate(samples = exp(samples)) %>% 
  
  ggplot(aes(x = samples, y = factor(dist, level = level_order))) +

  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(aes(fill = stat(x < 1)),
               .width = c(0.5, 0.95)
  ) +
  scale_fill_manual(values = c("gray90", "#919C4D")) +
  labs(x = "\nOdds Ratio",
       y = "Underlying prior\n",
       title = "Pooled data from all subgroups"
  ) +
  scale_x_continuous(breaks = seq(from = xlim_inf, to = xlim_sup, ticks_spacing)) +
  scale_y_discrete(expand = c(0, 0.3)) +
  theme(
    axis.title = element_text(size = 14),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 11),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = "none",
    plot.title.position = 'plot',
    plot.margin = margin(20, 25, 10, 20)
  ) +
  coord_cartesian(xlim = c(xlim_inf, xlim_sup))

p1 + plot_annotation(title = "Posterior distributions on mortality",
                          theme = theme(plot.title = element_text(size = 20)))

Interval bars depict 50% and 95% credible intervals

# Summary table
table_posterior = pooled_summary_table_death(all_posteriors_death)

# Use kableExtra to have a nice table

table_posterior %>% 
  kbl(booktabs = T, align = 'c') %>% 
  kable_styling(bootstrap_options = "striped", full_width = F) %>% 
  column_spec(5, bold = T, color = "#919C4D") %>%
  add_header_above(c(" " = 2,
                     "Probability of Benefit" = 3,
                     "Probability of Harm" = 1)) %>% 
  footnote(general = "OR = Odds Ratio")

		Probability of Benefit			Probability of Harm
Prior	Median OR	Pr(< 0.8)	Pr(< 0.9)	Pr(< 1)	Pr(> 1.1)
Non-informative	0.837	0.26	0.86	1.00	0
Evidence-based	0.846	0.19	0.83	1.00	0
Skeptical	0.856	0.15	0.79	0.99	0
Optimistic	0.832	0.27	0.89	1.00	0
Pessimistic	0.879	0.07	0.64	0.98	0
Note:
OR = Odds Ratio

Non-informative prior

xlabel = "\nOdds Ratio"
xlim_inf = 0.6
xlim_sup = 1.1

p1 = pooled_risk_difference_plot(
  recovery_death, # Data object

  ## start using quotes
  "#7EBAC2", # fill color code
  xlabel, # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.1 # X axis label
) 

p2 = pooled_cumulative_risk_difference_plot(
  recovery_death, # Data object
  xlabel, # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.05, # Spacing between ticks in X axis
  "mortality" # Outcome (within quotes)
  )

p1 + p2 + plot_annotation(
  title = "Posterior distribution: RECOVERY trial + Non-informative prior",
  subtitle = "Pooled data from all subgroups"
)

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the odds ratio (OR) is less than or equal to the effect size on the X‐axis.

Evidence-based prior

xlabel = "\nLog-Odds Ratio"

pooled_data_prior_posterior_plot(
  multiple_priors_death %>% # Output from pooled_normal_approximation()
    filter(type == "evidence-based"), 
  
  # start using quotes
           "#364D55", # Color for the prior
           "#7EBAC2", # Color for RECOVERY
           "#A65041", # Color for the posterior
           xlabel, # X axis label
           "Posterior distribution: RECOVERY + Evidence-based prior", # Title
           "Pooled data from all subgroups", # Subtitle
           # stop using quotes
           -0.50, # X axis inferior limit
           0.50, # X axis superior limit
           0.25 # Spacing between ticks in X axis
  ) 

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

xlim_inf = 0.6
xlim_sup = 1.1

p1 = pooled_posterior_difference_plot(
  multiple_priors_death %>% # Output from normal_approximation()
    filter(type == "evidence-based"), 

  ## start using quotes
  "#A65041", # fill color code
  "Odds Ratio", # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.1 # Spacing between ticks in X axis
) 

p2 = pooled_posterior_cumulative_plot(
  multiple_priors_death %>% # Output from normal_approximation()
    filter(type == "evidence-based"), 
  "Odds Ratio", # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.05, # Spacing between ticks in X axis
  "mortality" # Outcome (within quotes)
  )

p1 + p2 + plot_annotation(
  title = "Posterior distribution: RECOVERY trial + Evidence-based prior",
  subtitle = "Pooled data from all subgroups"
)

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the odds ratio (OR) is less than or equal to the effect size on the X‐axis.

Skeptical prior

xlabel = "\nLog-Odds Ratio"

pooled_data_prior_posterior_plot(
  multiple_priors_death %>% # Output from pooled_normal_approximation()
    filter(type == "skeptical"), 
  
  # start using quotes
           "gray50", # Color for the prior
           "#7EBAC2", # Color for RECOVERY
           "#A65041", # Color for the posterior
           xlabel, # X axis label
           "Posterior distribution: RECOVERY + Skeptical prior", # Title
           "Pooled data from all subgroups", # Subtitle
           # stop using quotes
           -0.50, # X axis inferior limit
           0.50, # X axis superior limit
           0.25 # Spacing between ticks in X axis
  ) 

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

xlim_inf = 0.6
xlim_sup = 1.1

p1 = pooled_posterior_difference_plot(
  multiple_priors_death %>% # Output from normal_approximation()
    filter(type == "skeptical"), 

  ## start using quotes
  "gray50", # fill color code
  "Odds Ratio", # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.1 # Spacing between ticks in X axis
) 

p2 = pooled_posterior_cumulative_plot(
  multiple_priors_death %>% # Output from normal_approximation()
    filter(type == "skeptical"), 
  "Odds Ratio", # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.05, # Spacing between ticks in X axis
  "mortality" # Outcome (within quotes)
  )

p1 + p2 + plot_annotation(
  title = "Posterior distribution: RECOVERY trial + Skeptical prior",
  subtitle = "Pooled data from all subgroups"
)

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the odds ratio (OR) is less than or equal to the effect size on the X‐axis.

Optimistic prior

xlabel = "\nLog-Odds Ratio"

pooled_data_prior_posterior_plot(
  multiple_priors_death %>% # Output from pooled_normal_approximation()
    filter(type == "optimistic"), 
  
  # start using quotes
           "#5CA881", # Color for the prior
           "#7EBAC2", # Color for RECOVERY
           "#A65041", # Color for the posterior
           xlabel, # X axis label
           "Posterior distribution: RECOVERY + Optimistic prior", # Title
           "Pooled data from all subgroups", # Subtitle
           # stop using quotes
           -0.50, # X axis inferior limit
           0.50, # X axis superior limit
           0.25 # Spacing between ticks in X axis
  ) 

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

xlim_inf = 0.6
xlim_sup = 1.1

p1 = pooled_posterior_difference_plot(
  multiple_priors_death %>% # Output from normal_approximation()
    filter(type == "optimistic"), 

  ## start using quotes
  "#5CA881", # fill color code
  "Odds Ratio", # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.1 # Spacing between ticks in X axis
) 

p2 = pooled_posterior_cumulative_plot(
  multiple_priors_death %>% # Output from normal_approximation()
    filter(type == "optimistic"), 
  "Odds Ratio", # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.05, # Spacing between ticks in X axis
  "mortality" # Outcome (within quotes)
  )

p1 + p2 + plot_annotation(
  title = "Posterior distribution: RECOVERY trial + Optimistic prior",
  subtitle = "Pooled data from all subgroups"
)

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the odds ratio (OR) is less than or equal to the effect size on the X‐axis.

Pessimistic prior

xlabel = "\nLog-Odds Ratio"

pooled_data_prior_posterior_plot(
  multiple_priors_death %>% # Output from pooled_normal_approximation()
    filter(type == "pessimistic"), 
  
  # start using quotes
           "#523C84", # Color for the prior
           "#7EBAC2", # Color for RECOVERY
           "#A65041", # Color for the posterior
           xlabel, # X axis label
           "Posterior distribution: RECOVERY + Pessimistic prior", # Title
           "Pooled data from all subgroups", # Subtitle
           # stop using quotes
           -0.50, # X axis inferior limit
           0.50, # X axis superior limit
           0.25 # Spacing between ticks in X axis
  ) 

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

xlim_inf = 0.6
xlim_sup = 1.1

p1 = pooled_posterior_difference_plot(
  multiple_priors_death %>% # Output from normal_approximation()
    filter(type == "pessimistic"), 

  ## start using quotes
  "#523C84", # fill color code
  "Odds Ratio", # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.1 # Spacing between ticks in X axis
) 

p2 = pooled_posterior_cumulative_plot(
  multiple_priors_death %>% # Output from normal_approximation()
    filter(type == "pessimistic"), 
  "Odds Ratio", # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.05, # Spacing between ticks in X axis
  "mortality" # Outcome (within quotes)
  )

p1 + p2 + plot_annotation(
  title = "Posterior distribution: RECOVERY trial + Pessimistic prior",
  subtitle = "Pooled data from all subgroups"
)

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the odds ratio (OR) is less than or equal to the effect size on the X‐axis.

Hospital discharge outcome

Now, let’s pool all the data and analyze the hospital discharge outcome. I believe odds ratio is the most appropriate effect size measure since there is great heterogeneity in baseline risk in each subgroup.

Prior

not_recovery %>% 
  filter(outcome == "discharge") %>% 
  select(-outcome) %>% 
  flextable() %>% 
  autofit()

Version 1: Random-Effect Meta-Analysis

prior_discharge =
  rma(
    # Outcome
    measure = "OR", # Odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = not_recovery %>% filter(outcome == "discharge"),

    slab = paste(trial_subgroup),
    method = "REML"
  )

# Visualize it
forest(prior_discharge,
       showweights = T,
       atransf=exp)

### Using escalc()

# Same results as rma() with random-effects

# effect_sizes =
#   escalc(
#     # Outcome
#     measure = "OR", # risk difference
#     
#     # Tocilizumab group
#     ai = events_toci,
#     n1i = total_toci,
# 
#     # Control group
#     ci = events_control,
#     n2i = total_control,
#     
#     data = not_recovery %>% filter(outcome == "discharge"),
# 
#     slab = paste(trial_subgroup)
#   )
# 
# res <- rma(yi, vi, data=effect_sizes) 
# 
# forest(res,
#        atransf=exp)

Version 2: Summing up all events + Fixed-Effect MA

total_discharge = not_recovery %>% 
  filter(outcome == "discharge") %>% 
  # Sum all events
  mutate(across(events_toci:total_control, ~ sum(.x))) %>% 
  # As all rows are now the same, pick 1
  slice(1) %>% 
  # Rename trial as "Total"
  mutate("trial_subgroup" = "Total")

prior_discharge_v2 =
  rma(
    # Outcome
    measure = "OR", # Odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = total_discharge,

    slab = paste(trial_subgroup),
    method = "REML"
  )

# Visualize it
forest(prior_discharge_v2,
       atransf=exp)

RECOVERY

recovery %>% 
  filter(outcome == "discharge") %>% 
  select(-outcome) %>% 
  flextable() %>% 
  autofit()

Version 1: Random-Effect Meta-Analysis

recovery_discharge =
  rma(
    # Outcome
    measure = "OR", # Odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = recovery %>% filter(outcome == "discharge"),

    slab = paste(trial_subgroup),
    method = "REML"
  )

# Visualize it
forest(recovery_discharge,
       showweights = T, 
       atransf=exp)

Version 2: Summing up all events + Fixed-Effect MA

total_discharge = recovery %>% 
  filter(outcome == "discharge") %>% 
  # Sum all events
  mutate(across(events_toci:total_control, ~ sum(.x))) %>% 
  # As all rows are now the same, pick 1
  slice(1) %>% 
  # Rename trial as "Total"
  mutate("trial_subgroup" = "Total")

recovery_total_discharge =
  rma(
    # Outcome
    measure = "OR", # Odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = total_discharge,

    slab = paste(trial_subgroup),
    method = "REML"
  )

# Visualize it
forest(recovery_total_discharge,
       atransf=exp)

Posterior distribution

Data overview

Regarding the hospital discharge outcome, odds ratio greater than 1 means benefit.

I am going to use data from Version 1 analyses. The effect size here will be log-odds ratio.

Regarding the skeptical prior, its mean is centered in 0, which is equal to log(1).

Regarding the optimistic prior, its mean is centered in 0.223, which is equal to log(1.25). The RECOVERY trial was designed to detect a relative risk reduction of 20% in the primary outcome (mortality). Thus, for the hospital discharge outcome, I decided to choose 0.223 for the mean since this is the exact opposite of -0.223, which in turn is equal to log(0.8). This latter value reflects the expected relative risk reduction mentioned above.

Regarding the pessimistic prior, its mean is centered in -0.223, which is equal to log(0.8). I mirrored the optimistic prior’s mean in the log-odds ratio scale.

multiple_priors_discharge %>% flextable() %>% autofit()

all_posteriors_discharge = pooled_tibble_all_posteriors(
  
  # Output from rma() with only RECOVERY
  recovery_discharge, 
  
  # Output from normal_approximation_multiple_priors()
  multiple_priors_discharge 
)

Summary

# X axis
xlim_inf = 1
xlim_sup = 1.6
ticks_spacing = 0.1

# Assure subgroup order
# https://stackoverflow.com/questions/12774210/
# how-do-you-specifically-order-ggplot2-x-axis-instead-of-alphabetical-order

level_order = c("Pessimistic",
                "Optimistic",
                "Skeptical",
                "Evidence-based",
                "Non-informative")

p1 = all_posteriors_discharge %>%
  
  # make it tidy for ggplot
  pivot_longer(
    cols = c("Non-informative":"Pessimistic"),
    names_to = "dist", values_to = "samples"
  ) %>%
  
  # exp() to transform log-odds ratio into Odds Ratio
  
  mutate(samples = exp(samples)) %>% 
  
  ggplot(aes(x = samples, y = factor(dist, level = level_order))) +

  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(aes(fill = stat(x > 1)),
               .width = c(0.5, 0.95)
  ) +
  scale_fill_manual(values = c("gray90", "#5891BF")) +
  labs(x = "\nOdds Ratio",
       y = "Underlying prior\n",
       title = "Pooled data from all subgroups"
  ) +
  scale_x_continuous(breaks = seq(from = xlim_inf, to = xlim_sup, ticks_spacing)) +
  scale_y_discrete(expand = c(0, 0.3)) +
  theme(
    axis.title = element_text(size = 14),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 11),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = "none",
    plot.title.position = 'plot',
    plot.margin = margin(20, 25, 10, 20)
  ) +
  coord_cartesian(xlim = c(xlim_inf, xlim_sup))

p1 + plot_annotation(title = "Posterior distributions on hospital discharge",
                          theme = theme(plot.title = element_text(size = 20)))

Interval bars depict 50% and 95% credible intervals

# Summary table
table_posterior = pooled_summary_table_discharge(all_posteriors_discharge)

# Use kableExtra to have a nice table

table_posterior %>% 
  kbl(booktabs = T, align = 'c') %>% 
  kable_styling(bootstrap_options = "striped", full_width = F) %>% 
  column_spec(3, bold = T, color = "#5891BF") %>%
  add_header_above(c(" " = 2,
                     "Probability of Benefit" = 3)) %>% 
  footnote(general = "OR = Odds Ratio")

		Probability of Benefit
Prior	Median OR	Pr(> 1)	Pr(> 1.2)	Pr(> 1.4)
Non-informative	1.346	1	0.96	0.28
Evidence-based	1.336	1	0.97	0.21
Skeptical	1.244	1	0.73	0.02
Optimistic	1.321	1	0.95	0.16
Pessimistic	1.172	1	0.34	0.00
Note:
OR = Odds Ratio

Non-informative prior

xlabel = "\nOdds Ratio"
xlim_inf = 1
xlim_sup = 1.6

p1 = pooled_risk_difference_plot(
  recovery_discharge, # Data object

  ## start using quotes
  "#7EBAC2", # fill color code
  xlabel, # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.1 # X axis label
) 

p2 = pooled_cumulative_risk_difference_plot(
  recovery_discharge, # Data object
  xlabel, # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.1, # Spacing between ticks in X axis
  "discharge" # Outcome (within quotes)
  )

p1 + p2 + plot_annotation(
  title = "Posterior distribution: RECOVERY trial + Non-informative prior",
  subtitle = "Pooled data from all subgroups"
)

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the odds ratio (OR) is less than or equal to the effect size on the X‐axis.

Evidence-based prior

xlabel = "\nLog-Odds Ratio"

pooled_data_prior_posterior_plot(
  multiple_priors_discharge %>% # Output from pooled_normal_approximation()
    filter(type == "evidence-based"), 
  
  # start using quotes
           "#364D55", # Color for the prior
           "#7EBAC2", # Color for RECOVERY
           "#A65041", # Color for the posterior
           xlabel, # X axis label
           "Posterior distribution: RECOVERY + Evidence-based prior", # Title
           "Pooled data from all subgroups", # Subtitle
           # stop using quotes
           -0.50, # X axis inferior limit
           0.50, # X axis superior limit
           0.25 # Spacing between ticks in X axis
  ) 

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

xlim_inf = 1
xlim_sup = 1.6

p1 = pooled_posterior_difference_plot(
  multiple_priors_discharge %>% # Output from normal_approximation()
    filter(type == "evidence-based"), 

  ## start using quotes
  "#A65041", # fill color code
  "Odds Ratio", # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.1 # Spacing between ticks in X axis
) 

p2 = pooled_posterior_cumulative_plot(
  multiple_priors_discharge %>% # Output from normal_approximation()
    filter(type == "evidence-based"), 
  "Odds Ratio", # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.1, # Spacing between ticks in X axis
  "discharge" # Outcome (within quotes)
  )

p1 + p2 + plot_annotation(
  title = "Posterior distribution: RECOVERY trial + Evidence-based prior",
  subtitle = "Pooled data from all subgroups"
)

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the odds ratio (OR) is greater than or equal to the effect size on the X‐axis.

Skeptical prior

xlabel = "\nLog-Odds Ratio"

pooled_data_prior_posterior_plot(
  multiple_priors_discharge %>% # Output from pooled_normal_approximation()
    filter(type == "skeptical"), 
  
  # start using quotes
           "gray50", # Color for the prior
           "#7EBAC2", # Color for RECOVERY
           "#A65041", # Color for the posterior
           xlabel, # X axis label
           "Posterior distribution: RECOVERY + Skeptical prior", # Title
           "Pooled data from all subgroups", # Subtitle
           # stop using quotes
           -0.50, # X axis inferior limit
           0.50, # X axis superior limit
           0.25 # Spacing between ticks in X axis
  ) 

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

xlim_inf = 1
xlim_sup = 1.6

p1 = pooled_posterior_difference_plot(
  multiple_priors_discharge %>% # Output from normal_approximation()
    filter(type == "skeptical"), 

  ## start using quotes
  "gray50", # fill color code
  "Odds Ratio", # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.1 # Spacing between ticks in X axis
) 

p2 = pooled_posterior_cumulative_plot(
  multiple_priors_discharge %>% # Output from normal_approximation()
    filter(type == "skeptical"), 
  "Odds Ratio", # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.1, # Spacing between ticks in X axis
  "discharge" # Outcome (within quotes)
  )

p1 + p2 + plot_annotation(
  title = "Posterior distribution: RECOVERY trial + Skeptical prior",
  subtitle = "Pooled data from all subgroups"
)

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the odds ratio (OR) is greater than or equal to the effect size on the X‐axis.

Optimistic prior

xlabel = "\nLog-Odds Ratio"

pooled_data_prior_posterior_plot(
  multiple_priors_discharge %>% # Output from pooled_normal_approximation()
    filter(type == "optimistic"), 
  
  # start using quotes
           "#5CA881", # Color for the prior
           "#7EBAC2", # Color for RECOVERY
           "#A65041", # Color for the posterior
           xlabel, # X axis label
           "Posterior distribution: RECOVERY + Optimistic prior", # Title
           "Pooled data from all subgroups", # Subtitle
           # stop using quotes
           -0.50, # X axis inferior limit
           0.50, # X axis superior limit
           0.25 # Spacing between ticks in X axis
  ) 

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

xlim_inf = 1
xlim_sup = 1.6

p1 = pooled_posterior_difference_plot(
  multiple_priors_discharge %>% # Output from normal_approximation()
    filter(type == "optimistic"), 

  ## start using quotes
  "#5CA881", # fill color code
  "Odds Ratio", # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.1 # Spacing between ticks in X axis
) 

p2 = pooled_posterior_cumulative_plot(
  multiple_priors_discharge %>% # Output from normal_approximation()
    filter(type == "optimistic"), 
  "Odds Ratio", # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.1, # Spacing between ticks in X axis
  "discharge" # Outcome (within quotes)
  )

p1 + p2 + plot_annotation(
  title = "Posterior distribution: RECOVERY trial + Optimistic prior",
  subtitle = "Pooled data from all subgroups"
)

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the odds ratio (OR) is greater than or equal to the effect size on the X‐axis.

Pessimistic prior

xlabel = "\nLog-Odds Ratio"

pooled_data_prior_posterior_plot(
  multiple_priors_discharge %>% # Output from pooled_normal_approximation()
    filter(type == "pessimistic"), 
  
  # start using quotes
           "#523C84", # Color for the prior
           "#7EBAC2", # Color for RECOVERY
           "#A65041", # Color for the posterior
           xlabel, # X axis label
           "Posterior distribution: RECOVERY + Pessimistic prior", # Title
           "Pooled data from all subgroups", # Subtitle
           # stop using quotes
           -0.50, # X axis inferior limit
           0.50, # X axis superior limit
           0.25 # Spacing between ticks in X axis
  ) 

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

xlim_inf = 1
xlim_sup = 1.6

p1 = pooled_posterior_difference_plot(
  multiple_priors_discharge %>% # Output from normal_approximation()
    filter(type == "pessimistic"), 

  ## start using quotes
  "#523C84", # fill color code
  "Odds Ratio", # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.1 # Spacing between ticks in X axis
) 

p2 = pooled_posterior_cumulative_plot(
  multiple_priors_discharge %>% # Output from normal_approximation()
    filter(type == "pessimistic"), 
  "Odds Ratio", # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  0.1, # Spacing between ticks in X axis
  "discharge" # Outcome (within quotes)
  )

p1 + p2 + plot_annotation(
  title = "Posterior distribution: RECOVERY trial + Pessimistic prior",
  subtitle = "Pooled data from all subgroups"
)

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the odds ratio (OR) is greater than or equal to the effect size on the X‐axis.

sessionInfo()

## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.3.4 patchwork_1.1.1  ggdist_2.4.0     metafor_2.4-0   
##  [5] Matrix_1.3-2     flextable_0.6.5  janitor_2.1.0    here_1.0.1      
##  [9] forcats_0.5.1    stringr_1.4.0    dplyr_1.0.5      purrr_0.3.4     
## [13] readr_1.4.0      tidyr_1.1.3      tibble_3.1.1     ggplot2_3.3.3   
## [17] tidyverse_1.3.1 
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-152          fs_1.5.0              lubridate_1.7.10     
##  [4] webshot_0.5.2         httr_1.4.2            rprojroot_2.0.2      
##  [7] tools_4.0.4           backports_1.2.1       bslib_0.2.4          
## [10] utf8_1.2.1            R6_2.5.0              DBI_1.1.1            
## [13] colorspace_2.0-1      withr_2.4.2           tidyselect_1.1.1     
## [16] curl_4.3.1            compiler_4.0.4        cli_2.5.0            
## [19] rvest_1.0.0           xml2_1.3.2            officer_0.3.18       
## [22] sass_0.3.1            scales_1.1.1          systemfonts_1.0.1    
## [25] digest_0.6.27         rmarkdown_2.7         svglite_2.0.0        
## [28] base64enc_0.1-3       pkgconfig_2.0.3       htmltools_0.5.1.1    
## [31] highr_0.9             dbplyr_2.1.1          fastmap_1.1.0        
## [34] rlang_0.4.11          readxl_1.3.1          rstudioapi_0.13      
## [37] jquerylib_0.1.4       generics_0.1.0        farver_2.1.0         
## [40] jsonlite_1.7.2        zip_2.1.1             distributional_0.2.2 
## [43] magrittr_2.0.1        Rcpp_1.0.6            munsell_0.5.0        
## [46] fansi_0.4.2           gdtools_0.2.3         lifecycle_1.0.0      
## [49] stringi_1.5.3         yaml_2.2.1            snakecase_0.11.0     
## [52] grid_4.0.4            crayon_1.4.1          lattice_0.20-41      
## [55] haven_2.4.1           hms_1.0.0             knitr_1.33           
## [58] pillar_1.6.0          uuid_0.1-4            reprex_2.0.0         
## [61] glue_1.4.2            evaluate_0.14         V8_3.4.0             
## [64] data.table_1.14.0     modelr_0.1.8          vctrs_0.3.8          
## [67] cellranger_1.1.0      gtable_0.3.0          assertthat_0.2.1     
## [70] cachem_1.0.4          xfun_0.22             broom_0.7.6          
## [73] viridisLite_0.4.0     memoise_2.0.0         ellipsis_0.3.2       
## [76] rdocsyntax_0.4.1.9000